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Monte  Carlo  simulations  are  being  used  to  investigate  the  structural  behavior  of 
bidisperse  colloidal  (silica)  liquids.  The  solids  volume  fraction  is  40%,  and  the 
bidispersity  is  characterized  by  two  parameters:  a,  the  particle  size  ratio,  and  p\  the 
relative  volume  fraction  of  small  particles  over  the  total  volume  fraction  of  solids.  Four 
different  size  ratios  are  used  (two,  three,  four,  and  five-to-one).  The  largest  particle  size  is 
0.6  um,  while  the  smallest  is  0.12  um.  The  variables  (electrolyte  concentration,  surface 
potential,  etc.)  of  the  potential  are  set  to  ensure  the  system  is  stable  and  non-flocculating. 

The  radial  distribution  function  g(r)  for  all  combinations  (small-small,  small- 
large,  and  large-large)  of  particles  along  with  the  osmotic  compressibility  is  calculated  for 
the  various  size  ratios  at  different  p  values.  There  are  two  interesting  changes  between  the 
monodisperse  and  bidisperse  suspensions.  First,  g(r)  for  the  large-large  particles  shows  an 
interesting  structural  change  in  that  the  HCP  (hexagonal  close  packing)  structure  for  the 
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monodisperse  large  particles  evolves  to  a  BCC-like  (body  centered  cubic)  structure  as  (3 
is  increased.  There  is  also  an  indication  of  phase  separation  which  may  be  a  consequence 
of  depletion  flocculation.  These  changes  appear  to  occur  for  all  size  ratios  that  have  been 
investigated.  Also,  the  HCP  (hexagonal  close  packing)  structure  for  the  monodisperse 
small  particles  is  completely  broken  up  when  added  to  the  large  particle  system.  Second, 
the  structural  change  appears  approximately  at  the  P  value  corresponding  to  the  minimum 
in  the  osmotic  compressibility.  The  structure  factor  is  also  calculated  from  the  Fourier 
transform  of  the  radial  distribution  function.  It  is  suggested  as  future  work  to  obtain  an 
effective  potential  between  large  particle  potentials,  treating  the  small  particle-laden  fluid 
as  a  continuum.  This  effective  potential,  which  can  be  obtained  using  an  inversion 
scheme,  may  indicate  an  attractive  force  between  large  particles  if  depletion  flocculation 
is  occurring. 
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CHAPTER  1 
INTRODUCTION 

Colloidal  suspensions  contain  microscopic  particles  suspended  in  a  liquid.  The 
size  of  these  particles  can  range  from  nanometers  to  a  few  micrometers.  The  particles  can 
be  biological  (viruses,  proteins,  and  bacteria),  as  well  as  manmade,  such  as  those  in 
paints,  cosmetics,  inks,  and  detergents.  Naturally  occurring  substances  such  as  clay, 
blood,  and  milk  can  also  be  considered  suspensions.  In  short,  understanding  colloidal 
suspensions  is  essential  in  many  areas. 

The  rheological  behavior  of  colloidal  suspensions  is  very  important  in 
determining  the  cost  effectiveness  of  processing  and  transporting  these  materials.  The 
higher  the  volume  fraction  of  solids,  the  higher  the  efficiency,  which  is  very  important  to 
industry.  However,  highly  concentrated  suspensions  can  be  very  viscous.  Therein  lies  the 
key  -  the  highest  throughput  while  maintaining  a  flowable  suspension. 

The  first  step  in  understanding  how  these  suspensions  behave  is  understanding 
how  the  particles  in  the  fluid  interact  with  one  another.  The  net  force  between  colloids  at 
equilibrium  consists  of  a  steric  repulsion  (if  coated  by  polymer),  electrostatic  repulsion, 
and  Van  der  Waals  attraction.  These  forces  are  dependent  on  the  physicochemical 
properties  of  the  fluid  medium  (pH,  electrolyte  concentration,  polymer  concentration, 
etc.)  as  well  as  the  particles  themselves.  For  instance,  at  high  electrolyte  concentrations, 
silica  particles  tend  to  agglomerate  due  to  the  attraction  caused  by  the  suppression  of  the 


1 


2 

electrical  double  layer.  When  electrolytes  are  present  in  low  concentrations,  the  systems 
tend  to  be  stable  because  of  electrostatic  repulsion.  The  latter  will  be  the  focus  of  this 
dissertation. 

The  microstructure  of  the  suspended  colloids  determines  the  fluidity  of  their 
suspensions.  Microstructure  is  determined  by  many  factors,  such  as  interparticle 
interaction,  polydispersity,  and  total  solids  volume  fraction.  Colloidal  suspensions  are 
found  in  gas-like,  liquid,  solid,  and  even  glassy  states.  Each  state  is  characterized  by 
distinct  microstructure  patterns. 

The  focus  of  this  project  is  to  examine  bidisperse  suspensions  consisting  of  two 
particle  sizes  as  a  starting  point  to  the  more  complicated  polydisperse  suspensions.  This 
work  focuses  on  systems  at  equilibrium  where  the  equilibrium  microstructure  is 
investigated  by  Monte  Carlo  simulations.  The  configurations  so  produced  can  be  used  for 
the  initialization  of  Stokesian  (Bossis  &  Brady  1984,  Brady  &  Bossis  1985, 1988) 
dynamic  simulations.  The  dynamic  simulations  study  directly  the  rheological  behavior  of 
the  suspension.  This  work  is  part  of  a  larger  effort  conducted  by  the  Engineering 
Research  Center  for  "Particle  Science  and  Technology"  at  the  University  of  Florida.  The 
effort  includes  experimental  studies  of  the  rheological  behavior  of  colloidal  suspensions 
and  also  focuses  on  bidisperse  suspensions. 

The  rest  of  this  thesis  has  been  organized  as  follows.  The  second  chapter  contains 
an  overview  of  prior  work  on  the  microstructure  of  colloidal  suspensions.  The  third 
chapter  focuses  on  the  microstructure  of  bidisperse  suspensions  as  revealed  by 
equilibrium  Monte  Carlo  simulations.  The  fourth  chapter  describes  our  efforts  to  deduce 
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an  effective  potential  between  large  particles  in  bidisperse  suspensions.  The  fifth  chapter 
discusses  our  results  of  dynamic  simulations  for  a  simple  shear  flow. 


CHAPTER  2 
LITERATURE  REVIEW 

2.1  Colloids  as  Model  Systems  for  Atoms 

In  the  description  of  equilibrium  properties  of  colloidal  suspensions  it  is  adequate 
to  adopt  a  coarse-grained  view  point  and  consider  the  liquid  medium  as  a  continuum  (van 
Megen  &  Snook  1984a).  The  fluid  affects  suspension  features  through  its  macroscopic 
properties  like,  the  dielectric  constant  and  density.  The  'potential  of  mean  force'(Poon  & 
Pusey  1995)  is  calculated  by  assuming  pair  additivity  for  the  pair  potentials,  which 
consist  of  a  Van  der  Waals  attraction  and  an  electrostatic  repulsion.  In  this  context  the 
equilibrium  thermodynamic  properties  for  the  colloidal  suspension  are  essentially  the 
same  as  those  of  an  atomic  system  with  an  equivalent  inter-atomic  potential.  Under 
certain  circumstances,  the  colloidal  particles  exhibit  the  same  structural  features  as  those 
found  in  atomic  systems  -  gas,  liquid,  crystal,  and  even  glasses.  A  colloidal  fluid  lacks  the 
long-range  order  found  in  a  colloidal  crystal,  but  is  characterized  by  well-defined  cells  of 
nearest  neighbors.  The  availability  of  well-characterized  models  for  colloidal  systems  has 
lead  to  a  growing  interest  in  using  colloids  as  model  systems  for  atoms. 

An  advantage  in  using  colloids  as  model  systems  for  atoms  is  the  capability  of 
controlling  the  interparticle  forces.  One  component  of  these  forces  is  typically  a  screened 
Coulomb  repulsion  of  the  form  (Verwey  &  Overbeek  1948) 
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Here  «p  is  the  particle  concentration,  each  particle  having  a  surface  charge  Zpe  and  a 
radius  ap.  The  dielectric  constant,  e,  is  a  macroscopic  property  of  the  fluid  medium;  kB  is 
the  Boltzmann  constant,  and  T  is  the  temperature.  The  number  density  na  and  charge  zae 
are  for  additional  ions  which  also  contribute  to  the  screening.  The  potential  given  in 
equation  [2-1]  is  called  the  Derjaguin-Landau-Verwey-Overbeek  (DLVO)  potential.  By 
varying  such  factors  as  the  electrolyte  concentration,  the  interparticle  forces  can  be 
manipulated.  The  force  as  well  as  its  range  can  be  controlled  by  changing  experimental 
conditions,  which  is  unfeasible  for  atomic  systems. 

Another  advantage  of  studying  colloids  over  atomic  systems  is  their  larger  size. 
The  size  of  some  particles  allows  for  direct  observation  by  video  microscopy.  Light 
scattering  techniques  also  allow  for  the  probing  of  particle  arrangements. 

In  normal  atomic  systems,  the  temperature  and  pressure  are  varied  to  study  phase 
transitions.  This,  however,  is  not  necessary  in  colloidal  suspensions.  The  ionic  strength 
plays  a  role  similar  to  that  of  temperature  in  atomic  systems.  As  the  ionic  strength 
increases  the  screening  length  decreases  causing  the  system  to  melt.  The  osmotic  pressure 
of  the  system  is  controlled  by  the  volume  fraction  of  solids. 
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2.2  Comparison  with  Atomic  Systems 

Colloidal  particles  in  suspension  are  at  equilibrium  with  the  fluid  medium  of 
specific  chemical  potential.  By  treating  the  colloids  as  "supramolecules"  and  the  fluid  as 
a  continuum,  the  statistical  mechanics  of  the  suspension  can  be  dealt  with  using  the  same 
methods  as  those  developed  for  atomic  liquids.  The  potential  between  the  particles 
corresponds  to  the  "bare"  potential  between  the  molecules  (McMillan  &  Mayerl945,  Vrij 
etal.  1978). 

The  length  scales  between  the  atomic  and  colloidal  systems  are  notably  disparate. 
This  makes  the  particle  density  differ  by  a  factor  of  ~  1010  (atomic  systems,  np  ~  1022  cm"3, 
whereas  colloids  np  ~  1012  cm"3).  The  elastic  constants  of  the  atomic  and  colloidal  systems 
differ  by  factors  of  the  same  order  (i.e.,  1010).  This  difference  can  be  rationalized 
qualitatively  by  the  fact  that  G  ~  Unp,  where  U  is  the  interaction  energy  between  the 
particles  (Crandall  &  Williams  1977).  In  the  two  systems,  the  energy  scale  U  is  in  the 
same  range  (~  10  eV),  as  can  be  expected,  since  colloidal  crystals  are  as  stable  as  the 
atomic  solids  at  room  temperature.  The  colloidal  solids  average  range  of  elastic  constants 
is  0.5  to  1000  dynes/cm2,  because  of  this,  it  is  very  easy  for  the  system  to  shear  melt.  This 
can  be  accomplished  by  a  simply  bumping  or  shaking  the  crystal. 

Colloids  experience  Brownian  dynamics  due  to  the  presence  of  the  fluid  medium, 
whereas  conventional  solids  undergo  ballistic  dynamics.  The  thermodynamic  properties 
of  the  system  are  not  changed  by  this  difference,  but  the  dynamics  are  affected.  The 
backflow  of  the  liquid  in  colloidal  crystals  overdamps  the  longitudinal  sound  modes,  and 
only  transverse  phonons  (quanta  of  sound  modes)  of  small  wavevector  are  underdamped 


7 

(Hurd  et  al.  1982).  In  typical  crystals,  however,  both  longitudinal  and  transverse  sound 
modes  are  capable  of  being  propagated. 

Even  though  there  has  been  improvement  in  the  production  of  monodisperse 
colloidal  particles,  most  colloids  are  still  mass  produced  with  a  distribution  in  the  size  of 
the  particles.  There  also  is  the  charge  polydispersity  to  consider.  While  both  are  common 
in  colloids,  there  is  nothing  similar  in  atomic  systems.  The  structure  and  dynamics  of  the 
particles  are  very  dependent  on  both  polydispersities  (size  and  charge)  (Dickinson  et  al. 
1981,  Dickinson  &  Parker  1985,  Barrat  &  Hansen  1986,  Pusey  1987). 

2.3  Theory  Behind  the  Distribution  Function 

The  correct  thermodynamics  of  liquids  cannot  be  predicted  by  the  virial  equation 
of  state  used  for  imperfect  gases.  Liquids  are  so  dense  that  many-body  interactions 
dominate.  Therefore,  different  approaches,  focused  on  the  prediction  of  the  radial 
distribution  function  g(r)  were  developed.  The  radial  distribution  function  can  be 
considered  a  component  that  multiplies  the  bulk  density  p  to  give  a  local  density  p(r)  = 
pg(r)  about  some  stationary  molecule.  There  are  two  key  reasons  why  g(r)  is  the  central 
quantity  of  liquid-state  theories.  First,  all  thermodynamic  functions  for  the  system  can  be 
expressed  as  integrals  over  g(r)  if  the  total  potential  energy  of  the  N-body  system  is  pair- 
wise  additive  meaning 

UN(r„...,rN)  =  £u(ry) 

i<J  [2-3] 

where  the  summation  is  over  all  pairs  of  molecules.  Second,  the  radial  distribution 

function  can  be  measured  experimentally  by  X-ray  diffraction. 
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For  any  liquid-state  theory  to  be  complete  an  equation  for  g(r)  is  needed.  First, 
h(r,2)  is  defined  as  g(r12)  - 1  and  is  a  measure  of  the  total  influence  of  molecule  1  on 
molecule  2  at  a  distance  r,2.  In  1914  Ornstein  and  Zernike  (reproduced  by  Frisch  & 
Lebowitz  1964)  proposed  the  separation  of  h(r12)  into  two  parts,  a  direct  part  and  an 
indirect  part.  The  direct  part  is  called  the  direct  correlation  function  c(r12).  The  indirect 
part  comes  from  the  contribution  of  all  other  molecules  in  the  system.  The  following 
equation  can  then  be  written 

h(r12)  =  c(r)2)  +  pjc(r13)h(r23)dr3 

The  above  expression  is  called  the  Ornstein-Zernike  equation  and  can  be  thought  of  as  the 
defining  equation  for  the  direct  correlation  function.  If  c(r)  were  written  in  terms  of  h(r) 
or  g(r)  then  equation  [2-4]  could  be  solved.  This  is  accomplished  by  introducing  an 
approximation  for  c(r).  c(r)  represents  the  direct  correlation  between  two  particles  and 
therefore  can  be  approximated  by  the  following  expression. 

c(r)  =  g,otai(r)  "  gfadtoaO")  [2-5] 
Before  proceeding  a  new  quantity  must  be  defined.  This  quantity  is  called  the  potential  of 
mean  force  w  and  is  the  interaction  between  two  molecules  a  fixed  distance  apart  while 
also  accounting  for  the  effect  of  other  molecules.  At  low  densities  the  potential  of  mean 
force  is  very  similar  with  the  intermolecular  potential  «(r).  The  following  expression 
serves  as  the  defining  equation  for  w 

««  -  <=*  [2-6] 
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P  is  equal  to  1/kT.  g,ota](r)  is  the  radial  distribution  function  g(r)  and  gindirect(r)  is  g(r) 
without  the  direct  interaction  «(r)  accounted  for;  i.e.  gindirect  =  exp[-P{w(r)-u(r)}].  Making 
these  substitutions  equation  [2-5]  becomes 

c(r;  -  e        -  e  [2-7] 
Two  more  functions  are  introduced  to  complete  the  formulation  of  the  Percus-Yevick 
equation;  the  Mayer  function7(r)  and  the  cavity  correlation  function  y(r).  J{t)  and  y(r)  are 
defined  below: 


/(r)  =  e-p"(r)  -  1 


y(r)  =  ep"(r)g(r) 


[2-8] 
[2-9] 


Rewriting  equation  [2-7]  in  terms  of  these  new  functions  gives 


c(r)  =  /(r)y(r)         (PY)  [2.1Q] 

The  result  of  substituting  the  equation  labeled  PY  into  the  Ornstein-Zernike  equation 
(equation  [2-4])  is  the  Percus-Yevick  (PY)  equation: 

y(r12)  =  1  +  pj/(r13)y(r13)h(r23)dr3 

To  restate,  the  Percus-Yevick  equation  is  obtained  by  assuming  an  expression  for  c(r), 
called  a  closure  relationship  (equation  [2-10]),  which  in  turn  is  substituted  into  the 
Ornstein-Zernike  equation.  The  PY  equation  is  then  solved  for  w(r).  This  allows  the 
determination  of  g(r)  and  all  thermodynamic  properties,  as  discussed  earlier.  We  must 
stress  that  the  PY  equation  is  noi  exact,  since  the  closure  (equation  [2-7])  is  only  an 
approximation.  "Statistical  Mechanics"  by  McQuarrie  (1973)  is  a  good  source  for  a  more 
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in  depth  discussion  of  the  Percus- Yevick  equation  and  integral  theories  of  liquids  in 
general. 

Although  the  Percus- Yevick  equation  was  first  derived  for  monatomic  liquids  it 
can  also  be  used  for  colloidal  suspensions  (most  work  has  been  done  for  hard  spheres). 
The  major  difference  being  the  use  of  an  interparticle  potential  instead  of  a  inter-atomic 
potential.  Experimentally,  light  scattering  (instead  of  X-ray  diffraction)  is  used  for 
colloidal  suspensions  because  of  the  size  of  the  particles. 

Compressibility  and  other  equilibrium  properties  for  hard  sphere  mixtures 
calculated  using  the  Percus- Yevick  integral  equation  have  been  compared  with  Monte 
Carlo  and  molecular  dynamics  data  by  Mansoori  et  al.  (1971).  They  found  there  was  good 
agreement  between  the  equation  of  state  and  the  machine  calculated  data.  The  PY 
equation  for  monodisperse  hard  spheres  was  solved  independently  by  Theile  (1963)  and 
Wertheim  (1963).  Lebowitz  (1964)  solved  it  for  the  case  of  a  mixture  of  hard  spheres  by 
the  method  of  functional  Taylor  expansion.  The  solution  is  in  good  accord  up  to  the  third 
power  in  density  with  the  virial  expansion  of  the  pressure  and  in  the  limiting 
monodisperse  case  reduces  to  that  previously  calculated  by  Thiele  and  Wertheim. 

Hard-sphere  perturbation  theory  has  been  compared  with  Monte  Carlo  simulations 
by  van  Megen  and  Snook  (1984b).  The  theoretical  results  are  obtained  by  dividing  the 
potential  energy  into  two  parts—potential  energy  of  an  unperturbed  (reference)  system  and 
that  of  the  perturbed  system-and  performing  a  perturbation  expansion  on  the  canonical 
partition  function  of  the  system  and  on  the  configurational  integral  in  particular  (for 
further  details  refer  to  McQuarrie).  First,  the  osmotic  compressibility  was  compared 
between  the  two  at  various  conditions.  The  conditions  were  manipulated  to  look  at  both 
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stable  systems  and  those  with  a  significant  secondary  minimum.  The  number  of  equal 
sized  particles  used  in  the  Monte  Carlo  simulations  were  32  or  108.  From  their  results, 
they  concluded  the  hard  sphere  perturbation  theory  can  be  used  to  accurately  predict 
thermodynamic  properties.  They  also  looked  at  how  well  the  radial  distribution  function 
calculated  using  perturbation  theory  compared  with  MC  results  and  found  little 
discrepancy  between  the  two.  Overall,  they  concluded  the  hard-sphere  perturbation  theory 
could  be  used  to  accurately  predict  the  structure  (g(r))  and  the  osmotic  pressure  of 
concentrated  colloidal  dispersions  whose  structure  is  essentially  determined  by  the  hard 
core  or  steeply  repulsive  part  of  the  interparticle  potential. 

The  radial  distribution  function  is  also  related  to  the  experimentally  measured 
structure  factor,  S(q),  by  (Allen  &  Tildesley  1994) 


where  q  =  (47in0/A.)sin9/2  is  the  magnitude  of  the  scattering  vector,  r^  is  the  refractive 
index  of  the  dispersion,  X  the  vacuum  wavelength  of  the  incident  radiation,  and  9  the 
scattering  angle.  The  low  q  region  of  S(q)  is  dominated  by  the  long  range  microstructure 
of  the  suspension  and  is  limited  by  the  instrumentation  whereas  the  large  q  region  is 
dominated  by  short  range  microstructure.  S(q)  illustrates  the  Fourier  components  of 
density  fluctuations  in  the  suspension. 

There  have  been  very  few  experimental  studies  which  focused  on  the 
microstructure  of  bidisperse  suspensions.  Bartlett  and  Ottewill  (1992)  have  performed 
light  scattering  on  concentrated  dispersions  of  sterically-stabilized  poly(methyl 
methacrylate)  particles  of  approximately  a  three-to-one  size  ratio  using  small-angle 
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neutron  scattering  experiments.  The  dispersing  medium  was  a  mixture  of  cis-decalin  and 
octane.  From  their  findings  they  concluded  the  possibility  of  small  particle  clustering  due 
to  an  attractive  force  caused  by  the  depletion  of  large  particles  from  in  between  the  small 
particles.  Duits  et  al.  (1991)  have  used  small-angle  neutron  scattering  experiments  to 
study  bidisperse  mixtures  of  different  size  silica  cores  coated  with  a  layer  of  octadecyl 
chains  (19.9  nm  and  52.7  nm).  The  dispersing  medium  was  cyclohexane.  The 
experiments  were  performed  at  three  total  volume  fractions,  0.22,  0.30,  and  0.40,  with 
equal  partial  volume  fractions.  Partial  structure  factors  were  obtained  using  the  small- 
angle  neutron  scattering  contrast  variation  method.  The  structure  factors  of  the 
monodisperse  suspensions  correlated  rather  well  with  those  calculated  from  the  hard- 
sphere  PY  approximations.  On  the  other  hand,  the  bidisperse  partial  structure  factors 
showed  a  possible  attraction  between  the  different  size  particles  which  is  thought  to  occur 
from  differences  in  the  grafting  density  on  both  particles.  The  main  contribution  from 
their  study  is  the  analysis  method  used  to  recover  the  partial  structure  factors.  Rodriguez 
and  Kaler  (1992a)  studied  the  equilibrium  structure  of  monodisperse  and  bidisperse  (141 
and  84  nm)  polystyrene  latex  particles  dispersed  in  bromoform  using  static  light 
scattering.  The  monodisperse  systems  structure  factor  was  predicted  with  good  agreement 
using  the  PY  hard-sphere  approximation.  The  bidisperse  partial  structure  factors  matched 
relatively  well  at  low  volume  fractions,  but  deviated  from  theory  as  the  volume  fraction 
increased.  The  results  also  deviated  as  the  volume  fraction  of  small  particles  was 
increased.  They  were  unable  to  give  any  reasons  for  the  difference. 

Although  there  has  been  considerable  work  done  on  bidisperse  hard  spheres,  few 
studies  have  used  simulations  to  investigate  the  structure  of  bidisperse  colloidal 


13 

suspensions.  The  main  contribution  of  our  work  has  been  the  relevance  of  structural 
changes  in  bidisperse  colloidal  suspensions  and  how  the  bidispersity  effects  other 
equilibrium  properties. 


CHAPTER  3 
EQUILIBRIUM  SIMULATIONS 

3.1  Introduction 

The  overall  goal  of  this  research  is  to  obtain  a  better  understanding  of  how 
polydispersity  affects  the  rheology  of  colloidal  suspensions  by  using  simulations  as  an 
investigative  tool.  Our  simulations  have  focused  on  bidisperse  suspensions  which  will 
serve  as  a  stepping  stone  to  the  polydisperse  case.  We  have  employed  Monte  Carlo 
simulations  first:  to  investigate  suspension  microstructure  under  equilibrium  conditions. 
The  next  step  will  be  to  study  suspension  microstructure  and  rheology  under  shear  flow 
conditions,  using  Stokesian  dynamics  (Bossis  &  Brady  1984,  Brady  &  Bossis  1985, 
1988).  In  this  chapter  we  present  our  findings  on  g(r)  of  concentrated  bidisperse 
suspensions. 

In  order  to  simulate  highly  concentrated  colloidal  suspensions,  yet  in  the  liquid 
regime,  the  total  volume  fraction  for  all  size  ratios  is  set  at  40%.  The  freezing  volume 
fraction  for  hard  spheres  is  49.4%  and  beneath  this  the  system  is  a  stable  fluid  (Hoover  & 
Ree  1968).  While  the  total  volume  fraction  of  particles  is  fixed,  the  composition  of  small 
and  large  particles  is  varied.  Following  the  well-established  Metropolis  Monte  Carlo 
method,  the  pair  correlation  functions  for  the  large-large  and  small-small  particles  are 
determined,  as  well  as  the  osmotic  compressibility  and  the  mean  square  displacement  of 
the  small  and  large  particles.  The  fluidity  of  these  suspensions  is  displayed  in  figure  3-1. 
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The  mean  square  displacement  shows  how  far  the  particles  have  moved  and  is  a  measure 
of  particle  diffusion.  The  pair  correlation  functions  also  display  the  liquid-like 
characteristics,  which  is  evident  by  the  weak  first  minimum  in  g(r)  for  the  monodisperse 
suspensions.  The  systems  begin  from  a  HCP  (hexagonal  close  packing)  configuration  to 
make  sure  the  suspensions  do  in  fact  melt  to  a  liquid.  Figure  3-2  reveals  that  this  is  in  fact 
the  case.  The  systems  are  also  started  from  random  initial  configurations. 

For  the  conditions  in  this  research,  the  monodisperse  suspensions  of  large  and 
small  particles  show  a  HCP  structure.  The  structure  mentioned  is  that  of  the  particles  in 
the  fluid  suspension.  However,  both  particles  (small  and  large)  in  the  bidisperse  mixture 
develop  an  arrangement  quite  different  from  that  of  the  monodisperse  case.  This  change 
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Figure  3-1  Mean-Square  Displacement  vs.  Monte  Carlo  simulation  cycles  for  large  (0.6 
um)  and  small  (0.2  um)  particles.  Both  of  which  are  monodisperse  suspensions.  (MSD 
has  been  scaled  by  the  square  of  the  respective  particle  radius.  The  coefficient  of 
correlation  (R2)  is  greater  than  0.99.) 


Figure  3-2  The  pair  correlation  function  for  two  monodisperse  cases.  The  open  squares 
are  for  the  0.2  jam  size  particles  and  the  closed  circles  are  for  the  0.6  um  size  particles. 
The  x-axis  is  scaled  by  the  diameter  of  the  particle  in  each  system. 


becomes  very  pronounced  when  p  is  about  15%.  p  is  the  relative  volume  fraction  of  small 
particles  over  the  total  volume  fraction  of  solids.  At  approximately  this  same 
composition,  the  compressibility  has  a  minimum. 

The  following  chapter  will  be  divided  into  four  sections.  The  first  section  will 
discuss  monodisperse  suspensions.  The  second  section  will  be  a  description  of  the 
simulation  method  where  a  brief  discussion  of  Metropolis  Monte  Carlo  will  be  provided. 
The  conditions  of  the  simulation  will  also  be  given.  The  following  section  will  cover  the 
results  on  g(r)  along  with  a  discussion,  and  the  final  section  will  be  a  summary. 
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3.2  Monodisperse  Suspensions 

In  order  to  put  our  findings  on  bidisperse  suspensions  in  context  we  performed 
preliminary  investigations  of  the  monodisperse  suspensions.  This  was  done  to  ensure  the 
systems  were  in  fact  in  the  liquid  regime.  As  mentioned  previously,  hard  sphere 
crystallization  begins  at  a  volume  fraction  of  49.4%  solids.  The  nominal  volume  fraction 
for  the  silica  particles  must  be  below  this  due  to  the  electrostatic  repulsion  of  the  potential 
between  the  particles.  This  repulsive  nature  effectively  increases  the  volume  fraction  of 
the  solids.  Two  forms  of  criteria  were  used  to  justify  the  fluidity  of  the  suspensions:  the 
radial  distribution  function  and  the  mean  square  displacement.  The  fact  that  the  first 
minimum  in  the  radial  distribution  function  is  far  from  zero  is  indicative  of  a  fluid  (figure 
3-2).  The  mean  square  displacement  displays  evidence  that  the  particles  diffuse  a 
substantial  distance  which  would  only  happen  in  a  liquid  (figure  3-1).  The  small  particles 
diffuse  less  due  to  the  increased  densification  caused  by  the  more  repulsive  nature  of  the 
small  particles.  For  these  reasons,  a  nominal  volume  fraction  of  40%  is  used.  At  this 
volume  fraction,  it  is  ensured  that  the  suspension  is  a  liquid  at  all  p  values  (bidisperse 
case).  It  also  should  be  mentioned  that  all  monodisperse  suspensions  are  arranged  in  a 
HCP  structure,  although  the  HCP  structure  is  less  pronounced  with  the  large  particles. 
This  is  discussed  further  in  section  3.4. 

3.3  Simulation  Method 

The  equilibrium  simulation  is  that  of  the  Metropolis  Monte  Carlo.  The  Monte 
Carlo  method  was  created  by  von  Neumann,  Ulam,  and  Metropolis  soon  after  the  Second 
World  War  to  look  at  the  diffusion  of  neutrons  in  fissionable  material.  'Monte  Carlo'  was 
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coined  by  Metropolis  on  account  of  the  broad  utilization  of  random  numbers  in  the 
computation.  The  initial  work  can  be  found  in  the  paper  by  Metropolis  and  Ulam  (1949). 
The  Metropolis  Monte  Carlo  was  initiated  by  Metropolis  et  al.  in  1953. 

The  Metropolis  Monte  Carlo  method  is  designed  to  produce  a  trajectory  in  phase 
space  which  samples  from  a  selected  statistical  ensemble.  The  examination  of  n  using  the 
Monte  Carlo  procedure  is  an  example  of  this  method  as  an  integration  tool.  The 
integration  is  accomplished  by  calculating  the  area  of  a  circle  with  unit  radius.  The  circle 
inscribed  in  a  square  is  shown  in  figure  3-3. 

Random  numbers  from  a  uniform  distribution  between  zero  and  one  are 
generated.  The  random  numbers  are  then  used  as  coordinates  for  what  is  called  a  shot. 
The  shot  is  considered  a  hit  if  its  distance  from  the  origin  (marked  O)  is  less  than  one. 
The  shaded  area  comprises  the  successful  hit  zone,  n  is  calculated  using  the  following 
formula 

4  x  Area  under  the  curve  CD       4  xhi, 

71  ~  =   — 

Area  of  the  square  ODBC         Tsho«  [3  jj 

The  number  of  hits  is  represented  by  xhit  and  the  total  number  of  shots  is  t^.  The 
exactness  of  the  estimate  is  dependent  on  the  number  shots  taken  and  is  of  the  order 
(Tshot)"l/2-  This  is  a  simple  example  of  the  usefulness  of  Monte  Carlo. 
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Simple  one-dimensional  integration  can  be  done  with  well  known  numerical 
methods,  but  Monte  Carlo  is  very  useful  for  multidimensional  integration  which  is 
commonly  seen  in  statistical  mechanics.  An  example  of  this  would  be  the  calculation  of 
the  configurational  integral 

Zwr  =  Jdrexp(-pV). 

[3-2] 

V  is  the  potential  energy  and  p  is  1/kT.  The  configurational  integral  is  the  potential 
contribution  to  the  partition  function.  The  method  is  computationally  inefficient  when 
states  are  sampled  from  a  uniform  distribution.  Metropolis  et  al.  (1953)  developed  a 
method  known  as  importance  sampling  to  overcome  this  limitation.  From  this  comes  the 
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coined  name  Metropolis  Monte  Carlo.  The  main  idea  behind  importance  sampling  is  to 
sample  preferentially  those  states  that  contribute  substantially  to  the  integral.  This  is 
accomplished  by  using  a  Markov  chain  of  states.  According  to  Allen  &  Tildesley  (1994), 
"a  Markov  chain  is  a  sequence  of  trials  that  satisfies  two  conditions: 

(a)  The  outcome  of  each  trial  belongs  to  a  finite  set  of  outcomes,  {T,,  T2, ... 
rm,rn, ...},  called  the  state  space. 

(b)  The  outcome  of  each  trial  depends  only  on  the  outcome  of  the  trial  that 
immediately  precedes  it." 

In  the  actual  simulation  this  is  achieved  by  moving  a  particle  a  random  distance  which  is 
not  to  exceed  a  specified  distance,  called  8rmax.  The  initial  state  is  m  and  the  new  state  is 
n.  The  energy  difference  (5V„J  is  calculated  between  the  two  states.  If  the  energy  is 
lower,  the  move  is  accepted,  if  the  energy  is  higher,  the  probability  exp(-p5V„J  is 
calculated.  A  random  number  chosen  from  a  uniform  distribution  between  zero  and  one  is 
compared  to  this  probability;  if  it  is  lower  the  move  is  accepted,  otherwise  the  state 
remains  unchanged. 

5rraax  can  be  chosen  in  a  variety  of  ways.  One  is  to  vary  6rmax  so  that  the  acceptance 
ratio  for  a  particle  to  move  remains  at  a  set  value.  Some  researchers  (Wood  &  Jacobson 
1959)  indicate  that  an  acceptance  ratio  of  0.1  maximizes  the  movement  through  phase 
space.  However,  most  researchers  agree  that  it  is  inefficient  to  reject  90%  of  the 
attempted  moves,  therefore,  an  acceptance  ratio  of  0.5  is  generally  used.  Alternatively, 
one  can  fix  5rmax  so  that  the  acceptance  ratio  is  close  to  the  desired  value.  The  latter  case 
with  an  acceptance  ratio  of  0.5  was  used  in  our  simulations. 
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Another  issue  is  choosing  which  particles  to  move.  The  particles  can  be  chosen  at 
random  or  sequentially  picked.  One  Monte  Carlo  cycle  consists  of  N  (number  of 
particles)  attempted  moves  whether  the  particles  are  selected  sequentially  (each  particle 
has  an  index  number)  or  randomly.  Picking  the  particles  sequentially  reduces  the  random 
number  generation  and  is  just  as  correct  in  sampling  the  states  (Hastings  1970).  Our 
simulation  picks  the  particles  sequentially. 

The  Metropolis  Monte  Carlo  method  is  applied  to  colloidal  suspensions  of  two 
different  size  particles  under  conditions  typical  of  aqueous  suspensions  of  silica  particles. 
The  choice  of  silica  particles  was  motivated  by  concurrent  experimental  efforts  in  our  lab. 
The  sizes  range  from  0.12  urn  to  0.6  urn.  A  canonical  ensemble  of  500  to  6912  particles 
is  used  along  with  periodic  boundary  conditions  for  the  simulation.  The  number  of 
particles  is  chosen  to  ensure  that  there  are  no  defects  in  the  system  -  meaning  if  the 
system  were  HCP  there  would  be  no  empty  spaces  for  particles  to  reside  without 
disrupting  the  structure.  The  formula  used  to  calculate  the  number  of  particles  is 

N  =  4  *  NC3.  p_3] 
N  being  the  number  of  particles  and  NC  any  positive  integer.  Most  simulations  start  from 
at  least  two  different  random  configurations  and  move  toward  an  equilibrium  state  by 
following  the  Metropolis  importance  sampling  criterion  (Kehr  &  Binder  1984,  Allen  & 
Tildesley  1994).  Some  simulations  are  also  initiated  from  a  HCP  structure  to  ensure  that 
the  suspension  is  fluid.  It  is  well  known  that  in  a  N-V-T  simulations,  spontaneous  melting 
of  a  superheated  solid  is  difficult  (Allen  &  Tildesley  1994).  Therefore,  the  observation  of 
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spontaneous  melting  guarantees  that  our  thermodynamic  conditions  place  the  system  well 
inside  the  liquid  region  of  its  phase  diagram. 

Pairwise  additivity  of  interparticle  forces  is  assumed  in  calculating  the  total 
energy  of  the  system.  It  is  also  assumed  that  the  interparticle  interactions  consist  of 
screened  electrostatic  repulsion  and  a  van  der  Waals  attraction  for  which  the  following 
expressions  are  used  (Mahanty  &  Ninham  1976): 


i  t  ZteW°b(M>a+M>b) 
Ur  - 


l  +  exp(-Kh)l 
xyl  +  vl    ll-exp(-Kh)J      1       HV  U 


[3-4] 
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CTa°b      ,      <v?b      hlj£  -(q«+qb) 

r2-(aa+crb)2    r2-(aa-ab)2    2  |r2-(aa-ab)2 


[3-5] 

Here  s0  is  the  permittivity  of  free  space  and  e  the  dielectric  constant  of  the  medium,  a  is 
the  particle  radius,  y  is  the  surface  potential,  and  h  is  the  minimum  surface-to-surface 
distance  between  two  particles.  The  subscripts  a  and  b  refer  to  the  small  and  large 
particles,  respectively,  k  is  the  reciprocal  of  the  Debye  length,  which  characterizes  the 
double-layer  thickness: 

k2  =  —  yv2 

8Ki  1  [3-6] 
Here  e  is  the  fundamental  charge  (1.6  x  1019  C),  z,  and  n,  are  the  valency  and  number 
concentration  of  ion  i,  and  k  is  the  Boltzmann  constant,  r  in  equation  [3-5]  is  the  center-to 
center  distance  between  particles  and  A  is  the  Hamaker  constant. 

While  equation  [3-4]  is  an  approximation  for  the  case  of  low  surface  potential 
(i.e.,  zevj/j  «  kT),  it  is  known  to  be  appropriate  when  y,  <  60  mV  and  Kaa  >  5  (Hogg  et 
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al.  1966).  The  total  potential  between  two  particles  is  simply  the  sum  of  the  electrostatic 
repulsion  and  the  van  der  Waals  attraction: 

uT  =  uR  +  uA  [3J] 

For  an  aqueous  suspension  of  silica  particles,  s  is  set  to  be  80,  and  the  Hamaker  constant 
A  is  8.5  x  10"21  J.  While  the  surface  potential  of  the  particles  may  vary  depending  on  the 
pH  value  of  the  suspension,  47  mV  is  used  for  both  small  and  large  particles.  It  is  further 
assumed  that  the  suspension  contains  1-1  electrolyte  (e.g.,  NaCl)  at  a  concentration  of  1 
mmol/1.  Under  these  conditions,  the  secondary  minimum  in  the  interparticle  potential  is 
negligible  and  the  particles  remain  stable  in  suspension.  In  addition,  the  Debye  length  (or 
the  thickness  of  the  diffuse  layer)  relative  to  the  particle  radius  is  smaller  than  1/13  for  the 
small  particles  (i.e.,  ko,  >  13).  It  is  even  smaller  for  the  larger  particles,  thus  justifying 
the  use  of  equation  [3-4]. 

Figure  3-4  represents  the  interparticle  potential  between  two  particles  calculated 
using  equations  [3-4]  -  [3-7].  Under  the  prescribed  conditions,  the  particles  are  repulsive 
with  the  primary  maximum  for  the  interparticle  potential  greater  than  200  kT  in  all  cases. 
While  the  Debye  length  is  short  compared  to  the  particle  radius,  it  is  long  enough  to 
suppress  the  secondary  minimum,  which  always  appears  due  to  long  range  van  der  Waals 
attraction.  This  is  displayed  in  figure  3-5.  The  lowest  minimum  is  found  between  large 
particles,  because  of  the  relative  Debye  length.  The  minimum  is  -0.43  kT,  which  is  too 
weak  to  stabilize  agglomerates.  Therefore,  the  suspensions  simulated  are  well-dispersed 
and  stable. 
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Since  the  particle  volume  fraction  is  high,  special  attention  should  be  given  in 
creating  initial  configurations  to  avoid  particle  overlap.  Two  different  approaches  are 
used  for  the  initial  configurations.  In  the  first  approach,  random  coordinates  are  assigned 
to  the  particles.  A  very  low  solid  fraction  is  assumed  initially  to  avoid  particle  overlap. 
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Figure  3-4  Interparticle  potential  in  kT  units  as  a  function  of  surface-to-surface  distance 
scaled  by  the  small  particle  diameter.  (0.2  urn  and  0.6  um  size  particles) 


As  the  simulation  proceeds,  the  solid  fraction  is  gradually  increased  by  enlarging  the 
particle  size  to  a  desired  value  (i.e.,  40%  volume  fraction),  thereby  obtaining  a  non- 
overlapping  configuration.  In  the  second  approach,  the  solid  volume  fraction  is  fixed  at 
the  outset  and  random  coordinates  are  assigned  to  the  particles  as  in  the  first  approach. 
Since  the  system  is  very  dense,  particles  do  overlap.  This  difficulty,  however,  can  be 
overcome  by  eliminating  the  primary  minimum  in  the  interparticle  potential  artificially. 
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Since  the  height  of  the  primary  maximum  is  higher  than  200  kT  (figure  3-4),  it  is  unlikely 
for  the  particles  to  overcome  the  barrier  to  form  a  cluster.  Thus,  the  physical  consequence 
of  the  simulation  is  not  altered  by  modifying  the  interparticle  potential  to  a  constant  at  the 
maximum  value  in  the  region  where  the  interparticle  distance  is  smaller  than  where  the 
primary  maximum  occurs.  This  elimination  of  the  primary  minimum  makes  the 
overlapping  initial  configurations  highly  unfavorable  and  leads  quickly  to  non- 
overlapping  configurations. 
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Figure  3-5  Details  of  figure  3-4  on  a  much  finer  energy  that  allows  viewing  of  the 
secondary  minimum. 


Starting  from  at  least  two  different  initial  configurations  obtained  by  either 
method,  all  simulations  continue  for  at  least  20,000  cycles.  Some  simulations  are  run  for 
400,000  cycles  to  ensure  that  the  systems  are  in  fact  in  equilibrium.  The  compressibility, 
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as  well  as  the  energy,  does  not  change  after  the  20,000  cycles.  The  pair  correlation 
function  calculated  after  20,000  cycles  is  also  the  same  as  any  obtained  after  that.  This  is 
evident  in  figure  3-6.  The  radial  distribution  function  after  the  first  20,000  cycles  is 
plotted  along  with  one  calculated  after  40,000  cycles  and  there  is  no  difference. 
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Figure  3-6  The  radial  distribution  function  plotted  versus  center-to-center  distance  scaled 
by  the  particles  diameter.  This  figure  is  for  two  simulations:  one  after  20,000  cycles  (open 
squares)  and  the  other  after  40,000  cycles  (closed  circles).  The  two  curves  are  practically 
on  top  of  each  other.  The  plot  is  for  a  0.6  um  monodisperse  particle  suspension. 

Typically,  the  osmotic  compressibility  (Z)  and  the  radial  distribution  function 
(g(r))  are  calculated  from  200  configurations  that  are  separated  by  100  cycles.  The 
configurations  used  are  sampled  after  the  suspension  has  reached  equilibrium.  The  mean 
square  displacement  (MSD)  of  the  small  and  large  particles  is  calculated  using 
configurations  at  intervals  of  50  cycles.  Equilibrium  properties  are  calculated  as  ensemble 
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averages  of  microscopic  quantities  according  to  the  standard  formulas  of  statistical 
thermodynamics  (Allen  &  Tildesley  1994).  For  example,  the  osmotic  compressibility  is 
(van  Megen  &  Snook  1975) 


z.  J*-  =  i-  3 


NkT  NkT 


(%)  where  <^N)  =  ^.V^lA 


[3-8] 

Here  N  is  the  total  number  of  particles  in  the  system,  r^  and  Uy  are  the  relative  position 
vector  from  particle  i  to  particle  j  and  the  pair  potential  given  by  equation  [3-7], 
respectively.  The  bracket  '<  >'  denotes  an  ensemble  average.  The  radial  distribution 
function  is  calculated  from 


«<r>  =  :=r(II*<'-«i)! 


norm 


[3-9] 

The  delta  function  represents  a  function  which  is  non-zero  at  set  separations.  This  creates 
a  histogram  which  is  normalized  by  norm.  Norm  is  a  function  of  the  volume  of  the 
section,  the  number  of  particles,  and  the  number  of  configurations  sampled. 

The  macroscopic  properties  of  colloidal  suspensions  are  determined  by  various 
physicochemical  properties  (electrolyte  concentration,  particle  size,  particle  volume 
fraction,  and  temperature,  etc.)  as  can  be  seen  by  equations  [3-4]  -  [3-9].  When  the  system 
is  bidisperse,  there  are  two  additional  important  variables  that  come  into  play.  One  is  the 
particle  size  ratio  (aa/ab)  and  the  other  is  the  relative  volume  fraction  of  small  particles  to 
the  total  particle  volume  fraction  (p).  The  affect  of  these  variables  will  be  discussed  in  the 
following  section. 
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3.4  Results  and  Discussion 

Figure  3-7  displays  the  radial  distribution  function  for  three  different  size  particles 
in  a  monodisperse  suspension.  The  graph  shows  that  all  monodisperse  suspensions 
arrange  themselves  in  a  HCP  structure.  This  is  evident  by  the  broad  secondary  peak.  As 
the  size  of  the  particles  is  decreased,  the  Debye  length  relative  to  the  particle  size 
increases.  This  increases  the  effective  solids  volume  fraction  of  the  suspension.  As  a 
result  of  this  effective  densification  a  stronger  HCP  structure  develops  for  the  0.12  and 
0.2  urn,  which  is  apparent  by  the  split  in  the  second  peak.  The  radial  distribution 
functions  for  the  other  particle  sizes  (0.15  urn,  0.30  urn)  used  in  our  simulations  are  very 
similar  to  those  of  the  0.12  and  0.2  um  particles,  respectively. 
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Figure  3-7  Radial  distribution  function  for  monodisperse  suspensions.  The  solids  volume 
fraction  is  40%. 
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As  p  (relative  volume  fraction  of  small  particles  divided  by  the  total  solids 
volume  fraction)  is  increased,  an  interesting  change  appears  in  the  pair  correlation 
function  for  both  size  particles.  First  of  all,  the  large-large  particle  radial  distribution 
function  starts  to  develop  a  peak  in  between  the  primary  and  secondary  maximums.  This 
is  evident  for  all  size  ratios  and  can  be  seen  in  figures  3-9  to  3-12.  The  new  peak  starts  to 
emerge  at  a  p  very  close  to  the  minimum  in  the  compressibility  and  continues  to  become 
more  pronounced  as  p  increases. 

The  structural  change  depicted  in  figures  3-9  through  3-12  may  be  due  to  a 
clustering  of  large  particles  which  results  in  phase  separation.  The  two  phases  being  large 
particle  rich  or  large  particle  lean.  The  relatively  high  magnitude  of  the  first  peak  (greater 
than  five  in  each  case)  in  the  large  particles  radial  distribution  function  in  the  bidisperse 
suspensions  is  indicative  of  large  particle  clustering.  The  clustering  of  large  particles  may 
be  explained  in  the  context  of  depletion  flocculation.  In  case  of  hard  spheres,  the 
depletion  flocculation  may  be  due  to  increased  entropy  when  large  particles  cluster. 
Because  the  particles  are  not  able  to  overlap,  there  are  regions  where  the  small  particles 
center  of  mass  cannot  reside.  This  zone  is  the  distance  of  the  small  particles  radius  around 
the  large  particles.  Overlapping  of  the  depletion  of  two  large  particles  provides  more  free 
volume  for  the  small  particles  in  the  suspension,  which  increases  the  entropy.  It  may  be 
also  viewed  as  depletion  flocculation  due  to  the  exclusion  of  small  particles  in  the 
depletion  zone,  causing  a  net  osmotic  force  pushing  the  large  particles  together.  The 
movement  of  the  second  peak  in  the  large  particles  radial  distribution  function  as  p  is 
changed  could  be  due  to  the  relative  size  of  the  Debye  layer.  As  the  small  particles 
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become  smaller  they  become  more  repulsive,  imposing  a  stronger  force  on  the  clusters  of 
large  particles,  causing  the  large  particles  to  come  closer  together.  Therefore,  as  a,  the 
size  ratio  is  increased  the  secondary  peak  in  the  large  particles  radial  distribution  function 
(figures  3-9  through  3-12)  shifts  to  the  left. 

This  structural  change  may  be  also  explained  by  the  formation  of  "triplets"  (figure 
3-8)  with  one  small  particle  between  two  large  particles  in  a  linear  "chain"  formation.  For 
a  size  ratio  of  five-to-one,  the  peak  should  be  at  approximately  1.2.  For  hard  spheres  it 
would  indeed  be  very  close  to  this,  but  for  colloids,  the  double  layer  effectively  adds  to 
each  particle  radius.  Therefore,  the  location  of  the  peak  will  be  at  an  increased  distance. 


Figure  3-8  Schematic  of  a  "triplet". 

As  more  small  particles  are  added  to  the  system,  so  arises  the  possibility  of  the  "triplet" 
formation  occurring  more  frequently.  This  leads  to  more  developed  secondary  peaks  at 
higher  p  values  (figure  3-13).  The  figure  is  for  a  system  where  the  size  ratio  is  five-to- 
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one,  but  the  increase  in  the  peak  height  for  higher  (3  values  is  common  to  all  size  ratios 
simulated. 

The  small-small  particle  radial  distribution  function  in  the  bidisperse  suspension 
(figure  3-14)  shows  no  apparent  structure,  indicating  a  high  fluidity  for  the  small 
particles.  It  is  interesting  to  note  that  the  HCP  structure  originally  shown  by  the 
monodisperse  suspension  is  completely  broken  once  the  small  particles  are  added  with 
the  large  particles.  The  lack  of  second  and  third  nearest  neighbor  peaks  in  that  radial 
distribution  function,  also  suggests  the  absence  of  any  substantial  clusters  of  small 
particles  (at  least  for  the  p  values  depicted).  This  is  also  evident  in  figure  3-15  which  is 
the  total  radial  distribution  function  for  all  particle  interactions  (small-small,  small-large, 
and  large-large)  of  the  bidisperse  system  (0.2  and  0.6  urn).  The  first  and  second  peaks  are 
for  small-small  and  small-large  interactions,  respectively.  The  similar  magnitude  in  the 
first  two  peak  heights  may  indicate  it  is  just  as  likely  for  a  small  particle  to  be  next  to  a 
large  particle  as  it  is  another  small  particle.  This  same  behavior  is  seen  for  other  size 
ratios  as  well.  From  these  results,  it  can  be  inferred  that  increasing  the  relative  volume 
fraction  of  small  particles  induces  a  structural  rearrangement  of  the  large  particles,  while 
the  small  particles  microstructure  is  completely  unstructured.  This  picture  remains  valid 
at  relatively  small  p  values,  e.g.  around  p's,  which  correspond  to  the  compressibility 
minimum. 

If  the  structural  change  is  due  to  depletion  flocculation,  it  may  be  confirmed  by 
calculating  the  effective  potential  between  large  particles  which  may  show  an  attractive 
force.  The  effective  potential  can  be  obtained  from  the  inversion  of  the  large-large 


32 

particle  structure  factor  which  is  discussed  in  chapter  four.  On  the  other  hand,  if  the 
structural  change  is  due  to  the  "triplet"  formation,  it  may  be  confirmed  by  looking  at  the 
experimentally  obtained  structure  factor. 
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Figure  3-9  Radial  distribution  function  for  monodisperse  (0.6  (am,  p  =  0.0)  and  large- 
large  particles  in  a  bidisperse  system  where  the  size  ratio  is  five-to-one. 


33 


Figure  3-1 1  Radial  distribution  function  for  monodisperse  (0.6  um,  p  =  0.0)  and  large- 
large  particles  in  a  bidisperse  system  where  the  size  ratio  is  three-to-one. 
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Figure  3-13  Radial  distribution  function  for  large-large  particles  in  two  different  systems 
with  the  same  size  ratio  (five-to-one),  but  different  p  values. 
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Figure  3-15  Total  radial  distribution  function  for  bidisperse  system  (0.2  and  0.6  jam),  p 
equals  0.10.  g(r)  is  scaled  by  the  small  particles  diameter.  The  first,  second,  and  third 
peaks  represent  small-small,  small-large,  and  large-large  interactions,  respectively. 
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Figure  3-16  Osmotic  compressibility  (Z)  plotted  versus  p  for  various  a  values.  The  size 
of  the  particles  in  parenthesis  is  in  units  of  um. 

In  figure  3-16,  the  osmotic  compressibility  (Z),  calculated  from  our  data  using 
equation  [3-8],  is  shown  as  a  function  of  p  at  four  different  size  ratios.  Figure  3-17  shows 
the  same  thing  with  the  low  p  value  region  enlarged.  The  large  particle  radius  is  held 
constant  while  the  small  particles  radius  is  changed  accordingly.  The  compressibility 
decreases  as  p  increases,  reaches  a  minimum  in  each  case,  and  attains  a  maximum  at  p  = 
1.0  (i.e.,  monodisperse  small  particles).  The  position  of  each  minimum  is  shifted  to  the 
left  (lower  p)  as  a  (the  size  ratio)  is  increased.  The  minimum  also  deepens  as  a  increases. 
The  large  value  of  Z  at  p  =  1 .0  is  due  to  the  fact  that  the  relative  double  layer  thickness  to 
the  particle  size  is  much  larger  with  the  small  particles  than  with  the  large  particles. 
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It  should  be  pointed  out  that  this  non-linear  behavior  of  Z  with  P  was  predicted 
previously  for  bidisperse  mixtures  of  hard  spheres  by  Mansoori  et  al.  (1971).  In  figure  3- 
1 8,  Z  vs.  p  plot  is  given  for  various  values  of  particle  size  ratio  (a  =  oVcrJ.  This  plot  was 
obtained  using  the  analytic  expression  of  Mansoori  et  al.  (1971),  which  is  based  on 
Lebowitz's  solution  to  the  Percus-Yevick  equation  for  the  pair  correlation  function 
(Mansoori  et  al.  1971,  Lebowitz  1964).  The  analytic  result  has  been  shown  to  be 
consistent  with  the  Monte  Carlo  simulation  results  of  Rotenberg  (1965)  and  Smith  and 
Lea  (1963).  The  same  trend  has  also  been  observed  experimentally  using  polystyrene 
latex  particles  in  bromoform  (Rodriguez  &  Kaler  1992a). 
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Figure  3-17  Blow  up  of  low  p  values  from  figure  3-16. 
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Figure  3-18  Osmotic  compressibility  (Z)  of  bidisperse  hard  sphere  mixture  at  a  total 
solids  volume  fraction  of  40%  (plot  of  equation  [7]  in  Mansoori  et  al.  1971). 

Using  an  effective  particle  radius,  it  is  possible  to  compare  the  osmotic 
compressibility  calculated  from  Mansoori 's  analytical  expression  and  the  osmotic 
compressibility  obtained  from  the  Monte  Carlo  simulations.  The  effective  particle  radius 
may  be  calculated  in  two  different  ways.  The  first  and  simplest  way  is  to  consider  it  as  the 
position  of  the  first  maximum  in  the  radial  distribution  function  of  the  monodisperse 
suspension.  The  second  way  is  to  estimate  it  by  the  following  integral: 

HSD  =  R(MAX)  +  r  [l-exp(-U/kT)]dr. 

JR(MAX)  J  [\0] 

HSD  (hard  sphere  diameter)  is  the  effective  diameter.  R(MAX)  corresponds  to  the 
location  of  the  primary  maximum  in  the  potential  and  R(I)  is  the  approximate  location  the 
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potential  becomes  zero.  The  effective  radius  allows  for  the  calculation  of  an  effective 
volume  fraction.  The  relative  Debye  length  is  larger  for  small  particles,  causing  the 
effective  size  ratio  to  be  less  than  the  original  size  ratio.  The  compressibilities  obtained 
from  our  simulations  and  Mansoori's  analytical  expression  using  the  effective  volume 
fraction  are  compared  for  two  size  ratios  (a  =  3  &  5)  in  figures  3-19  and  3-20.  The 
effective  hard  sphere  radius  used  was  that  equivalent  to  the  maximum  position  of  the 
radial  distribution  function.  Apparently,  it  gave  the  better  results.  It  is  interesting  to  note 
that  the  location  where  the  minimum  in  the  osmotic  compressibility  occurs  matches  very 
closely  with  the  results  obtained  from  Mansoori's  analytical  expression.  However,  there 
is  some  discrepancy  in  the  magnitude  of  the  osmotic  compressibility. 
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Figure  3-19  The  osmotic  compressibility  obtained  from  simulation  and  Mansoori's 
analytical  expression  are  shown  versus  p.  The  original  volume  fraction  and  size  ratio  are 
40%  and  three-to-one,  respectively. 
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Figure  3-20  The  osmotic  compressibility  obtained  from  simulation  and  Mansoori's 
analytical  expression  are  shown  versus  p.  The  original  volume  fraction  and  size  ratio  are 
40%  and  five-to-one,  respectively. 

In  calculating  the  compressibility  of  the  bidisperse  mixture  of  hard  spheres,  only 
the  value  of  the  pair  correlation  function  at  the  point  of  contact  plays  a  role. 
Consequently,  the  implications  of  microstructural  changes,  such  as  the  new  features  in  the 
radial  distribution  function,  have  been  given  little  attention.  The  decrease  in  the  osmotic 
compressibility  seems  to  be  associated  with  the  structural  change  previously  described, 
suggesting  that  the  new  structure  results  in  lower  free-energy  at  the  given  volume 
fractions.  Therefore,  a  higher  particle  volume  fraction  may  be  attainable  before  a  solid 
state  is  reached.  In  fact,  it  has  been  shown  experimentally  that  a  higher  value  of 
maximum  packing  fraction  is  associated  with  a  smaller  value  of  osmotic  compressibility 
for  bidisperse  suspensions  (Rodriguez  &  Kaler  1992a).  This  structural  change  in 
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association  with  the  minimum  in  the  osmotic  compressibility  is  one  of  the  new  findings 
of  our  work. 

Figures  3-1  and  3-21  represent  the  mean  square  displacement  (MSD)  of  large  and 
small  particles  in  monodisperse  and  bidisperse  suspensions,  respectively.  Since  the  x-axis 
of  the  plot  is  not  time,  but  the  number  of  Monte  Carlo  simulation  cycles  (MCC),  the 
actual  value  of  particle  diffusivity  cannot  be  determined  from  these  results.  In  general,  the 
Monte  Carlo  method  may  not  be  reliable  for  the  quantitative  calculation  of  dynamic 
properties,  even  if  a  definite  time  interval  were  to  be  associated  with  each  Monte  Carlo 
cycle.  It  is  due  to  the  fact  that  the  hydrodynamic  interactions  are  completely  neglected  in 
the  equilibrium  Monte  Carlo  simulations.  Although,  they  are  irrelevant  for  the  calculation 
of  equilibrium  microstructure,  they  will  affect  particle  dynamics  even  in  a  suspension  at 
equilibrium.  Nevertheless,  the  Monte  Carlo  results  can  be  trusted  to  distinguish 
qualitatively,  between  liquid-like  and  solid-like  particle  mobility.  In  the  latter  case, 
particles  are  expected  to  move  only  a  small  fraction  of  their  radius  over  the  time  scale  of 
the  simulation.  The  linear  dependence  of  MSD  on  MCC  observed  in  the  simulations 
shows  clearly  that  the  suspension  simulations  are  in  a  liquid  state.  In  addition,  the  relative 
magnitude  of  the  slope  of  each  plot  should  represent  the  particle  diffusivity  in  a  relative 
sense.  It  is  apparent  from  figures  3-1  and  3-21  that  the  diffusivity  of  the  small  particles  in 
the  bidisperse  mixture  is  larger  by  two  orders  of  magnitude  than  that  in  the  monodisperse 
suspension,  despite  the  same  total  solid  fraction  of  40%.  The  large  particles  however, 
have  a  decrease  in  their  diffusivity  in  the  bidisperse  case.  This  may  be  explained  by  the 
fact  that  the  small-large  particle  interactions  are  more  repulsive  making  it  more  difficult 
for  the  large  particles  to  move  around.  The  lower  mobility  of  large  particles  can  also  be 
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realized  if  they  form  clusters.  Thus,  it  may  be  concluded  from  these  figures  that  the 
suspensions  are  not  glassy  but  liquid-like  at  the  total  particle  volume  fraction  of  40%  and 
that  the  mobility  of  the  small  particles  in  a  bidisperse  mixture  can  be  significantly  larger 
than  that  in  a  respective  monodisperse  suspension. 
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Figure  3-21  Mean-Square  Displacement  vs.  Monte  Carlo  simulation  cycle  for  large  and 
small  particles  in  a  bidisperse  suspension,  p  =  0.1 1  and  the  size  ratio  is  two-to-one.  (MSD 
has  been  scaled  by  the  square  of  the  respective  particle  radius.  The  coefficient  of 
correlation  (R2)  is  greater  than  0.99.) 
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Although  the  osmotic  compressibility  is  an  equilibrium  property,  it  provid 
qualitative  information  on  dynamic  properties,  such  as  zero  shear  rate  viscosity.  The 
linear  behavior  described  in  figures  3-16  and  3-17  appears  to  be  associated  with  the 
structural  changes  reflected  in  the  pair  correlation  function.  The  low  compressibility  and 
the  high  mobility  of  the  small  particles  in  the  bidisperse  suspension  suggests  that  the 
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shear  viscosity  of  the  bidisperse  mixture  may  be  lower  than  the  monodisperse 
suspensions  of  either  large  or  small  particles  at  the  same  value  of  total  particle  volume 
fraction.  In  fact,  such  a  trend  has  been  observed  experimentally  (Rodriguez  &  Kaler 
1992b). 

3.5  Summary 

It  has  been  shown  by  Monte  Carlo  simulations  that  the  equilibrium  structure  of 
bidisperse  colloidal  suspensions  can  be  quite  different  from  that  of  monodisperse 
suspensions.  While  the  stable  monodisperse  suspensions  in  the  present  investigation 
exhibited  a  HCP  structure  at  the  particle  volume  fraction  of  40%,  the  bidisperse  mixture 
showed  a  structural  change  as  the  relative  volume  fraction  of  the  small  and  large  particles 
was  varied.  The  structural  change  was  also  accompanied  by  a  decrease  in  the  osmotic 
compressibility.  In  addition,  the  mean-square  displacement  of  both  large  and  small 
particles  indicated  that  the  mobility  of  the  particles  could  be  substantially  higher  in  a 
bidisperse  suspension  than  in  a  monodisperse  suspension. 

The  results  suggest  that  the  maximum  packing  fraction  (i.e.,  the  maximum 
particle  volume  fraction  before  a  transition  to  a  glassy  or  crystalline  state)  of  bidisperse 
suspensions  can  be  higher  than  that  attainable  with  monodisperse  suspensions.  The 
investigation  also  indicates  that  the  shear  viscosity  of  the  bidisperse  mixture  may  be 
lower  than  that  of  the  monodisperse  suspensions  of  either  large  or  small  particles  at  the 
same  value  of  total  particle  volume  fraction.  Our  results  and  conclusions  are  consistent 
with  experimental  observations. 


CHAPTER  4 
STRUCTURE  FACTOR 

4.1  Introduction 

The  structure  factor  presented  in  Chapter  2  can  directly  probe  the  microstructure 
of  colloidal  suspensions  as  it  can  be  measured  experimentally  using  radiation  scattering 
techniques,  such  as  x-ray,  neutron,  and  light  scattering.  The  structure  factor  is  calculated 
from  the  Fourier  transform  of  the  radial  distribution  function.  For  low  q,  the  structure 
factor  can  directly  be  calculated  from  the  simulation  by  a  direct-averaging  procedure 
developed  by  Frenkel  et  al.  (1986).  In  theory,  the  radial  distribution  function  can  be 
calculated  from  the  experimentally  obtained  structure  factor  by  a  Fourier  transform,  but 
in  practice  this  transformation  procedure  is  limited  by  the  noise  in  the  low  q  region.  The 
inversion  of  the  structure  factor  can  also  provide  an  effective  potential  (Rajagopalan  & 
Rao  1997)  for  simulated  colloidal  suspensions.  The  derivation  of  an  effective  potential 
for  the  large-large  particle  interactions  in  bidisperse  suspensions  is  the  main  focus  of  this 
chapter.  This  effective  potential  views  the  small  particles  as  a  part  of  the  fluid  medium, 
and  represents  the  overall  influence  of  the  small  particles  contribution  on  the  large-large 
particle  interactions.  The  following  section  will  address  briefly  the  relation  between  the 
structure  factor  and  g(r),  as  well  as  technical  aspects  of  the  structure  factor  calculation 
from  simulation  data.  It  will  also  expand  on  the  inversion  scheme  developed  by 
Rajagopalan  and  Rao  (1997). 


44 


45 

4.2  Theory 
The  local  particle  density  is  defined  as 

p(r)  =  Jftr-r,), 

i=l 

where  the  Fourier  components  of  the  density  p(r)  are 

Pq  =  Jexp(-iq»r)p(r)dr 


[4-1  j 


[4-2] 


=  ZexP("i(l#ri) 

i-i 

The  static  structure  factor  (called  previously  structure  factor)  is  defined  as 

[4-3] 

The  structure  factor  S(q)  characterizes  the  Fourier  components  of  the  density  fluctuations 
in  the  suspension.  If  the  system  is  homogeneous,  the  structure  factor  can  be  drafted  as  the 
Fourier  transform  of  g(r)  (Hansen  and  McDonald  1986): 


1 


N  N 


j    /  N  N 

=  1  +  S  J  J  exp[iq  •  Ci-  - 1-  )]5(r  -  )5(r  -  r^drdr 
=  1  +  pjexp(-iq«r)g(r)dr 

Conversely,  g(r)  is  the  Fourier  transform  of  S(q): 

P8(r)  =  (^p- JexPCi*!  •  r)[SCq|)  -  l]dq 


[4-4] 


[4-5] 
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For  isotropic  systems,  S(q)  can  be  written  in  terms  of  the  wavenumber  q  =  |q|.  With  this 

S(q)  =  1  +  47tpfr2g(r)  ^-^dr. 

qr  [4-6] 

Rewriting  equation  [4]  using 

h(q)  =  fexp(-iq»r)h(r)dr 

J  [4-7] 

gives 

S(q)  =  1  +  (27i)3p8(q)  +  ph(q) 

where  h(r),  the  total  correlation  function,  is  equal  to  g(r)  -  1.  Neglecting  the  forward 
scattering  of  radiation  in  experimentation,  S(q)  can  be  written  as 

S(q)  =  1  +  pfexp(-iq»r)h(r)dr. 

J  [4-9] 

In  Monte  Carlo  simulations,  S(q)  may  be  calculated  in  two  different  ways.  One 
way  is  to  obtain  S(q)  from  the  Fourier  transform  of  the  radial  distribution  function  using 
equation  [4-6].  The  S(q)  from  this  method  tends  to  be  very  noisy  in  the  low-q  region  due 
to  the  finite  size  of  the  system.  This  can  be  overcome  by  using  the  second  method,  which 
is  a  direct  averaging  of  [exp(-iq»r)]  during  the  simulation  by  using  the  first  line  in 
equation  [4-4]  (Salgi  et  al.  1992). 

The  inversion  scheme  used  to  obtain  the  pair  potential  from  S(q)  will  be  discussed 
in  the  following  paragraphs.  For  a  more  detailed  explanation  refer  to  Rajagopalan  &  Rao 
(1997).  h(r)  in  equation  [4-9]  is  related  to  the  direct  correlation  function  c(r)  by  the 
Ornstein-Zernike  equation  which  was  described  in  Chapter  2  (equation  [2-4]).  This 
equation  can  be  solved  using  the  exact  closure  relation: 
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c(r)  =  -pu(r)  +  h(r)  -  ln[h(r)  +  l]  +  B(r),  [41Q] 

where  B(r)  is  the  bridge  function  and  p  =  l/kBT.  There  is  no  way  of  analytically  obtaining 
B(r),  so  approximations  are  used.  One  of  these  approximations,  the  Hypernetted-Chain 
(HNC)  closure,  is  the  result  of  setting  the  bridge  function  equal  to  zero.  The  previously 
discussed  Percus-Yevick  (PY)  closure  (equation  [2-10]  in  Chapter  2)  is  achieved  by 
linearizing  the  HNC  closure  with  respect  to  [h(r)  -  c(r)]. 

For  any  inversion  scheme,  the  general  idea  is  to  use  the  Ornstein-Zernike  (OZ) 
equation  along  with  equation  [4-10]  to  capture  u(r)  from  S(q)  for  a  known  density  p. 
Using  the  OZ  equation  along  with  the  general  closure  relation  (equation  [4-10])  it  is 
possible  to  relate  c(r)  to  S(q)  by 


c(q)  =  - 
P 


1-  1 


S(q) 


[4-11] 

If  there  was  an  analytical  expression  for  the  bridge  function  it  would  be  possible,  using 
the  above  expressions,  to  obtain  u(r)  from  S(q),  but  since  there  isn't,  an  iterative  method 
is  used  to  attain  u(r). 

The  iterative  method  begins  by  writing  the  closure  relation,  equation  [4-10], 
relative  to  an  unknown  potential  u'(r)  as  follows: 


pu(r)  =  In 


y'(r) 


— — —  +  [h(r)-h'(r)]  -  [c(r)-c'(r)]  +  A[u(r),u'(r)], 

'  [4-12] 


where  the  prime  represents  the  quantities  corresponding  to  u'(r),  and  y'(r)  is  the  cavity 
function  given  by 

y'(r)  =  g'(r)e<*«>.  ^ 


Acquiring  the  correct  potential  u(r)  is  now  a  matter  of  obtaining  u'(r)  such  that  A[u(r), 
u'(r)]  ->  0,  where  A[u(r),  u'(r)]  =  [B(r)  -  B'(r)].  The  unknown  potential  u'(r)  is  replaced 
with  a  hard  sphere  potential  ud(r),  where  the  hard  sphere  diameter  d  is  now  the  iterative 
parameter.  Using  the  hard  sphere  potential  allows  the  use  of  the  extremely  accurate 
parameterized  form  for  the  hard-sphere  cavity  function  yd(r)  (Henderson  &  Grundke 
1975).  The  initial  estimate  of  d  is  obtained  using  the  relation  between  the  isothermal 
compressibility  and  S(q  =  0). 

>(PP) 


S(0) 


dp 


[4-14] 

One  can  use  the  Carnahan-Starling  equation  to  relate  the  pressure  P  to  the  diameter  d. 
The  potential  in  accord  with  this  diameter  is  then  split  into  two  components:  a  core  part 
uc(r)  and  a  perturbation  u„(r)  using  the  Weeks-Chandler-Andersen  criterion  (Andersen  et 
al.  1972). 

u(D  =  u.(r)  ♦  „„(r)  ^ 
Then,  using  the  core  part  of  the  potential  along  with  a  criterion  developed  by  Lado  (1982) 

fM)[e-K(r)_e-K(r)]dr  =  Q 

8d  [4-16] 
an  updated  diameter  is  calculated.  The  new  diameter  is  put  to  use  in  equation  [4-12]  to 
correct  u(r).  This  process  is  continued  until  d  and  u(r)  converge. 

4.3  Results 

In  this  section  the  structure  factor  calculated  from  the  Fourier  transform  of  the 
radial  distribution  function  is  compared  with  the  structure  factor  recovered  directly  from 
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the  simulation  using  direct  averaging.  The  structure  factor,  S(q)  for  both  the  small-small 
and  large-large  particles  in  the  bidisperse  suspension  is  shown.  Also,  the  results  for  the 
potential  obtained  from  the  inversion  of  the  structure  factor  are  given.  The  systems 
investigated  are  the  same  as  those  previously  discussed  with  a  set  volume  fraction  of  40 
%  and  a  size  ratio  of  three-to-one  (0.2  &  0.6  urn). 

As  stated  previously,  the  direct  averaging  is  only  needed  in  the  low  q  region,  so 
the  following  plots  will  not  show  the  whole  q  region,  but  just  the  region  of  direct 
averaging.  Figure  4-1  is  the  structure  factor  calculated  from  both  methods  for  the  large 
particles  in  a  bidisperse  suspension.  This  particular  simulation  ran  for  40,000  cycles,  and 
the  structure  factor  is  calculated  using  the  last  20,000  cycles.  Figure  4-2  shows  the  same 
comparison  of  the  two  different  methods  of  calculation  of  the  structure  factor  for  a 
monodisperse  (0.6  um)  suspension. 
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Figure  4-1  Structure  factor  for  large-large  particles  in  bidisperse  suspension.  The  particle 
size  ratio  is  three-to-one  (0.2  &  0.6  um),  and  p  is  0.10. 
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Figure  4-2  Structure  factor  for  a  monodisperse  (0.6  urn)  suspension  calculated  from  the 
Fourier  transform  of  g(r)  and  from  direct  averaging. 


The  structure  factor  for  the  large-large  and  small-small  particles  in  the  bidisperse 
suspensions  are  given  in  figures  4-3  through  4-7.  The  structure  factors  in  these  figures  are 
calculated  from  the  Fourier  transform  of  g(r).  Figure  4-3  shows  S(q)  for  both 
monodisperse  suspensions.  It  is  interesting  to  note  the  increase  in  the  first  peak  height  for 
the  0.2  urn  monodisperse  suspension.  Figures  4-4  and  4-5  depict  the  structure  factor  for 
the  large-large  particles  in  the  bidisperse  suspension.  The  magnitude  of  the  first  peak 
decreases  as  p  increases.  Figures  4-6  and  4-7  display  the  structure  factor  for  the  small- 
small  particles  in  the  bidisperse  suspension  with  the  p  value  indicated  in  each  case. 
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Figure  4-3  Structure  factor  calculated  from  the  Fourier  transform  of  the  radial  distribution 
function  for  both  monodisperse  suspensions. 
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Figure  4-4  The  structure  factor  for  large-large  particles  in  a  bidisperse  suspension  for 
various  p  values. 
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Figure  4-5  Same  figure  as  4-4  showing  a  blow  up  of  the  low  q  region. 
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Figure  4-6  The  structure  factor  for  the  small-small  particles  for  the  bidisperse  suspension 
for  various  p  values. 
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Figure  4-7  Same  as  figure  4-6  with  the  low  q  region  blown  up. 


.  The  inversion  of  the  structure  factor  to  obtain  the  effective  potential  was  carried 
out  for  a  monodisperse  large  particle  suspension.  The  inversion  of  structure  factor  for  the 
monodisperse  suspension  should  reproduce  the  pair  potential  which  was  used  in  the 
simulation.  Thus,  the  procedure  with  the  monodisperse  system  should  serve  as  a  way  to 
check  the  accuracy  of  the  inversion  method.  In  figure  4-8,  the  pair  potential  for  a 
monodisperse  large  particle  suspension  calculated  using  the  inversion  method  is 
compared  with  the  original  one  used  for  the  simulation  to  obtain  g(r).  While  the  position 
of  the  secondary  minimum  apparently  matches  well,  a  significant  discrepancy  is  shown  in 
the  magnitude  of  the  pair  potential.  Despite  several  trials,  we  could  not  resolve  the 
discrepancy  at  this  point,  and  for  that  reason  we  did  not  apply  the  inversion  method  to  the 
bidisperse  system.  Although  this  attempt  was  not  successful,  this  method  will  be  useful  to 
confirm  whether  the  structural  changes  discussed  in  chapter  three  for  the  bidisperse 
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suspensions  are  associated  with  the  depletion  flocculation  or  the  "triplet"  formation,  once 
the  difference  in  the  inversion  method  is  resolved. 

9.00  ,  


8.00 


Figure  4-8  Comparison  of  the  potential  obtained  from  the  inversion  of  the  monodisperse 
structure  factor  (effective  potential)  and  the  potential  used  in  the  simulation  (input 
potential).  The  interparticle  potential  is  in  kT  units  and  interparticle  distance  is  scaled  by 
the  particle  diameter. 


CHAPTER  5 

DYNAMIC  SIMULATION  FOR  SHEAR  VISCOSITY 
5.1  Introduction 

Rheological  behavior  of  colloidal  suspensions  is  a  function  of  many  parameters. 
Some  of  those  being  particle  size,  polydispersity,  dispersing  medium,  interparticle 
interactions,  and  particle  concentration.  It  would  be  beneficial  to  have  some  correlation 
between  the  equilibrium  properties  discussed  previously  and  the  rheological  behavior  of 
the  suspension.  We  have  investigated  the  rheology  of  the  suspension  using  the  Stokesian 
dynamics  simulation  method  developed  by  Brady  (Brady  &  Bossis  1988,  Durlofsky  et  al. 
1987).  In  the  following  section,  the  Stokesian  dynamics  simulation  is  discussed  followed 
by  our  results  obtained  from  the  simulations. 

5.2  Stokesian  Dynamics 

This  method  focuses  on  a  finite  system  of  rigid  particles  where  the  particle 
Reynolds  number,  Uap/\i,  is  much  less  than  one  so  that  the  inertia  is  negligible.  t/is  the 
characteristic  velocity,  a  the  characteristic  particle  size,  p  the  fluid  density,  and  p.  the 
fluid  viscosity.  The  particles  are  suspended  in  an  unbounded  Newtonian  fluid  where  their 
motion  is  governed  by  the  Stokes  equations 

V.u  =  0, 
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The  method  begins  with  the  integral  solution  for  Stokes  flow  (Ladyzhenskaya  1963) 
which  gives  the  velocity  of  any  particle  or  at  any  point  in  the  fluid: 

"'(x)  =  "r(x) "  iirti  W-MfrW, 

u.  (x)  is  the  velocity  field  in  the  absence  of  the  particles,  Sa  is  the  surface  of  particle  a,  y 
is  the  position  vector  on  the  particle  surface,  and  x  is  a  field  point  in  the  fluid-particle 
continuum.  Jtj  is  the  free-space  Green's  function  and  is  often  called  the  Stokeslet  or  Oseen 
tensor,  and  can  be  written  as 

8  rr 

Mr)  -  JL  +  A 

r       r  [5-3] 
where  r  =  x  -  y,  r  =  \r\.f/y)  is  the  force  density  at  the  point  y  on  the  surface  of  the 
particle.  The  integration  as  shown  by  the  summation  is  to  be  done  over  all  N  particle 
surfaces.  The  force  density  can  be  written  in  terms  of  the  fluid  stress  tensor  a*  by 

fj(y)  =  <**(yK(y).  [5.4] 

nk(y)  is  the  surface  normal  vector  pointing  into  the  fluid. 

The  force  density  (in  equation  [5-2])  on  the  surface  of  each  particle  is  then 
expanded  in  a  series  of  moments  about  the  center,  x"  of  each  particle: 


U((x)  _  ur(x)  =  __L|[js3(x_x)//v)d5y 


°yk  -I  [5-5] 

The  zeroth  moment  of  the  force  density  is  the  total  force,  while  the  first  moment  has  both 
antisymmetric  and  symmetric  parts  (the  couplet  (or  torque)  L  and  stresslet  S, 
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respectively).  The  stresslets  can  be  thought  of  as  the  contribution  to  the  bulk  stress  due  to 
the  existence  of  the  particles.  Using  Faxen  formula  for  spheres  (Batchelor  &  Green  1972) 
and  truncating  equation  [5-5]  after  the  first  moment,  the  grand  mobility  matrix  M  can  be 
created,  which  correlates  the  velocity  (both  translational  and  angular)  and  the  rate  of 
strain  of  each  particle  with  the  force,  torque,  and  stresslet  of  all  TV  particles  in  the  system: 


[5-6] 


"u-u00" 

"F" 

=  M« 

-E" 

S_ 

M  can  further  be  divided  into  four  submatrices: 


M  = 


[5-7] 


oo 

U  -  U  represents  both  the  translational  and  angular  velocity,  -  E"  is  the  bulk  rate  of 
strain,  and  for  simple  shear  is  given  as 


E°°  =  — 


0  1  0 

1  0  0 
0  0  0 


[5-8] 

F  is  the  vector  representing  the  force  and  torque  exerted  by  the  particles  on  the  fluid,  and 
S  corresponds  to  the  particle  stresslets.  MUF  relates  particle  velocities  to  forces, 
velocities  to  stresslets,  MEF  rate  of  strain  to  forces,  and  MH  rate  of  strain  to  stresslets. 

The  inverse  of  M  is  the  resistance  matrix  R.  The  inversion  of  the  mobility  matrix 
incorporates  the  many-body  problem  at  the  level  of  stresslets  and  forces  as  pointed  out  by 
Brady  &  Bossis  (1988).  The  resistance  matrix  relates  the  forces,  torques,  and  stresslets 
exerted  by  the  particles  on  the  fluid  to  the  particle  velocities  and  the  rate  of  strain. 
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"F" 

~u-ir~ 

S 

-E" 

[5-9] 

Even  though  the  resistance  matrix  accounts  for  many-body  interactions,  it  does  not 
represent  the  lubrication  force  properly  unless  all  moments  are  included  in  the  mobility 
matrix.  Lubrication  forces  are  due  to  the  thin  layer  of  viscous  fluid  that  separates  nearly 
touching  particles.  These  forces,  due  to  their  short  range  nature,  are  basically  two-body 
interactions,  and  can  be  accounted  for  in  a  pairwise-additive  manner  in  the  resistance 
matrix  using  the  exact  expression  known  for  the  two  sphere  problem  (Arp  &  Mason 
1977,  Jeffrey  &  Onishi  1984,  Kim  &  Mifflin  1985).  This  addition  is  represented  by  R2B 
(for  two-body  resistance).  However,  since  the  far-field  interactions  have  already  been 
included  in  the  inversion  of  the  mobility  matrix  M,  it  must  be  subtracted  off  from  R2B  so 

as  not  to  count  them  twice.  The  two-body  interactions  that  are  to  be  subtracted  off  (R*B ) 
are  calculated  from  the  inversion  of  a  two-sphere  mobility  matrix  which  is  approximated 
to  the  same  level  as  M  (i.e.,  truncated  after  the  first  moment).  Then  the  resistance  matrix 
that  includes  both  the  near-field  lubrication  and  the  far-field  many  body  interactions  can 
be  presented  as 

R  =  M"1  +  R2B  -  R^  [5.10] 

In  the  calculation  of  the  viscosity,  the  bulk  stress  <E>  is  needed.  The  bulk  stress  is 
determined  by  an  average  over  the  volume  ^containing  the  N particles  and  is  presented 


below: 


(l)  =  IT  +  2r1E-+^{(sH)  +  (sp)  +  (sB)} 

[5-1 1] 
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The  relative  viscosity  is  simply  (Chang  &  Powell  1993) 


H  xy  [5-12] 
In  equation  [5-11],  IT  is  an  isotropic  term  that  is  irrelevant.  The  addition  of  particles  to 
the  Newtonian  fluid  adds  three  new  contributions  to  the  bulk  stress.  The  first  term,  (SH>, 
is  a  mechanical  or  contact  stress  transmitted  by  the  fluid  due  to  the  shear  flow  and  is  the 
stresslet  mentioned  previously.  (Sp>  is  an  "elastic"  stress  due  to  the  interparticle  forces 
and  the  last  term  <SB>  is  the  contribution  from  the  Brownian  motion.  The  contributions  are 
presented  below: 

(sh)  =  (m-  .[-b--m,.fJ)  [5_n] 

(Sp)  =  (xFp) 

W      V     1  [5-14] 
Equation  [5-13]  is  obtained  from  equations  [5-6]  and  [5-7]  by  solving  for  S  the  stresslets. 
Fp  in  equation  [5-14]  is  the  non-hydrodynamic  forces  acting  on  the  particles  (e.g., 
colloidal  forces). 

5.3  Results 

The  Stokesian  dynamics  simulation  has  been  applied  to  a  stable  bidisperse 
suspension  described  in  chapter  three  (i.e.,  silica  particle  suspension  with  47  mV  surface 
potential  and  1  mM  electrolyte  concentration).  Since  this  dynamics  simulation  is  very 
time  consuming,  it  has  been  applied  to  a  monolayer  system  in  which  the  particles  are 
restricted  to  a  single  plane  forming  a  monolayer  of  particles.  This  reduces  the  degree  of 
freedom  from  1 1  to  6  per  particle  and  provides  a  suitable  numerical  setting,  because  it 
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reduces  the  computation  time  while  preserving  the  necessary  physics  in  the  plane  of 
shear.  For  this  stable  system  it  is  not  necessary  to  include  the  lubrication  forces,  as  the 
particles  do  not  come  close  to  each  other  due  to  electrostatic  repulsion.  This  also  reduces 
the  computation  time  significantly. 

The  bidisperse  suspension  consists  of  spherical  particles  of  0.2  jam  and  0.6  jam  in 
radius.  Both  monodisperse  and  bidisperse  suspensions  have  been  investigated.  The 
viscosities  calculated  for  various  shear  rates  are  displayed  below  for  both  monodisperse 
and  bidisperse  cases.  The  viscosity  for  the  bidisperse  case  is  shown  to  be  lower  than  both 
monodisperse  cases.  This  decrease  in  viscosity  with  the  bidisperse  suspension  apparently 
corresponds  well  with  the  minimum  in  the  osmotic  compressibility  for  the  same 
bidisperse  system  (chapter  three).  The  shear  thinning  behavior  for  both  monodisperse 
systems  may  indicate  the  presence  of  equilibrium  microstructure  by  colloidal  forces  (i.e., 
electrostatic  repulsion  in  this  case).  As  the  shear  rate  is  increased  further,  the  equilibrium 
structure  is  broken  and  hydrodynamic  forces  become  dominant,  leading  to  a  high  shear 
rate  viscosity  limit.  The  Newtonian  behavior  shown  by  the  bidisperse  suspension 
suggests  that  there  is  no  significant  equilibrium  structure  formation  at  this  area  fraction  of 
40  %.  It  is  believed  at  higher  area  fractions  the  shear  thinning  behavior  would  be  shown 
by  the  bidisperse  case  as  well. 
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Figure  5-1  Viscosity  versus  shear  rate  for  0.2  urn  and  0.6um  monodisperse  suspensions 
and  a  bidisperse  mixture  of  the  two.  The  p  value  for  the  bidisperse  suspension  is  0.1  and 
the  area  fraction  is  40%  solids. 


Figure  5-2  shows  the  initial  configuration  for  the  0.6  urn  monodisperse 
suspension  and  two  'snapshots'  after  being  sheared  for  0.97  seconds  and  3.79  seconds  at 
a  shear  rate  of  10  s1.  The  instantaneous  viscosity  of  the  suspension  is  shown  in  figure  5-3. 
Since  the  shear  viscosity  depends  on  the  instantaneous  configuration  of  particles,  it 
fluctuates  as  shown  in  this  figure.  While  figure  5-2  shows  an  evolution  of  the  particle 
configuration  from  a  structured  configuration  at  t  =  0  seconds  to  a  less  structured  one  at  t 
=  3.79  seconds,  figure  5-3  suggests  that  the  evolution  of  structure  may  repeat  in  a  cyclic 
manner.  Various  'snapshots'  for  the  bidisperse  case  are  shown  in  figure  5-4,  and  the 
instantaneous  viscosity  versus  time  is  shown  in  figure  5-5.  While  it  may  be  difficult  to 
comment  on  the  structural  evolution  in  figure  5-4,  we  may  note  that  the  shear  viscosity  of 
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the  bidisperse  suspension  (figure  5-5)  is  smaller  than  that  of  the  monodisperse  one  (figure 
5-3)  at  all  times. 


Initial  t  =  0.97  seconds  t  =  3 .79  seconds 
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Figure  5-2  Individual  'snapshots'  of  the  system  (0.6  um,  monodisperse)  as  it  is  being 
sheared.  The  instantaneous  viscosity  versus  time  for  this  simulation  is  shown  in  figure 
5-3. 
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Figure  5-3  The  instantaneous  viscosity  versus  time  is  plotted  for  the  monodisperse 
suspension  of  0.6  urn  silica  particles.  The  area  fraction  is  40%  and  the  shear  rate  is  10s1. 
Time  is  scaled  by  the  shear  rate. 
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Figure  5-4  'Snapshots'  at  various  times  are  shovvri  for  the  bidisperse  suspension  (0.2  (am 
and  0.6  urn).  The  area  fraction  is  40%  solids  and  the  shear  rate  is  1  s"1. 
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Figure  5-5  Viscosity  versus  time  for  the  bidisperse  suspension  (0.2  urn  and  0.6  urn).  The 
shear  rate  is  1  s"1  and  the  area  fraction  is  40  %  solids.  Time  is  scaled  by  the  shear  rate. 


5.4  Summary 

The  results  of  Stokesian  dynamics  simulation  show  that  the  stable  bidisperse 
suspensions  do  in  fact  have  a  lower  viscosity  than  both  monodisperse  suspensions.  This 
was  anticipated  from  our  equilibrium  results  which  showed  a  minimum  in  the  osmotic 


compressibility  for  the  bidisperse  suspension.  Since  low  osmotic  compressibility  is  an 
indication  of  more  favorable  packing  of  particles,  the  bidisperse  suspension  may  be  able 
to  accommodate  higher  volume  fractions  of  solid  particles.  Thus,  at  the  same  solid 
fraction,  the  viscosity  of  the  bidisperse  suspension  may  be  lower  than  that  of  the 
monodisperse  suspension.  The  more  favorable  packing  of  particles  can  be  a  consequence 
of  the  structural  change  described  in  chapter  3. 


CHAPTER  6 
CONCLUSIONS 

Monte  Carlo  simulations  were  used  to  investigate  the  structural  behavior  of 
bidisperse  colloidal  (silica)  liquids.  The  solids  volume  fraction  was  40%,  and  the 
bidispersity  was  characterized  by  two  parameters:  a,  the  particle  size  ratio  and  p,  the 
relative  volume  fraction  of  small  particles  over  the  total  volume  fraction  of  solids.  Four 
different  size  ratios  were  used  (two,  three,  four,  and  five-to-one).  The  largest  particle  size 
was  0.6  urn,  while  the  smallest  was  0.12  um.  The  variables  (electrolyte  concentration, 
surface  potential,  etc.)  of  the  potential  were  set  to  ensure  the  system  was  stable  and  non- 
flocculating. 

The  radial  distribution  function  g(r)  for  all  combinations  (small-small,  small- 
large,  and  large-large)  of  particles  along  with  the  osmotic  compressibility  was  calculated 
for  the  various  size  ratios  at  different  P  values.  There  were  two  interesting  changes 
between  the  monodisperse  and  bidisperse  suspensions.  First,  g(r)  for  the  large-large 
particles  showed  an  interesting  structural  change  in  that  the  HCP  (hexagonal  close 
packing)  structure  for  the  monodisperse  large  particles  evolved  to  a  BCC-like  (body 
centered  cubic)  structure  as  p  was  increased.  There  is  also  an  indication  of  phase 
separation  which  may  be  a  consequence  of  depletion  flocculation.  These  changes 
appeared  for  all  size  ratios  that  were  investigated.  Also,  the  HCP  (hexagonal  close 
packing)  structure  for  the  monodisperse  small  particles  was  completely  broken  up  when 
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added  to  the  large  particle  system.  Second,  the  structural  change  appeared  approximately 
at  the  p  value  corresponding  to  the  minimum  in  the  osmotic  compressibility.  The 
structure  factor  was  also  calculated  from  the  Fourier  transform  of  the  radial  distribution 
function.  It  is  suggested  as  future  work  to  obtain  an  effective  potential  between  large 
particle  potentials,  treating  the  small  particle-laden  fluid  as  a  continuum.  This  effective 
potential,  which  can  be  obtained  using  an  inversion  scheme  developed  by  Rajagopalan 
and  Rao  (1997),  may  indicate  an  attractive  force  between  large  particles  if  depletion 
flocculation  is  occurring. 
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